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Eric B. Gregory, 1 ' 2 Christine T. H. Davies, 1 ^ Iain D. Kendall, 1 Jonna Koponen, 1 Kit Wong, 1 Eduardo 
Follana, 3 Elvira Gamiz, 4 G. Peter Lepage, 5 Eike H. Muller, 6 Heechang Na, 7 and Junko Shigemitsu 7 

(HPQCD collaboration) 

1 School of Physics and Astronomy, The Kelvin Building, University of Glasgow, Glasgow G12-8QQ, UK 
2 Department of Physics, University of Cyprus, PO Box 20357, 1678 Nicosia, Cyprus 
Departamento de Ftsica Teorica, Universidad de Zaragoza, 
CI. Pedro Cerbuna 12, E- 50009 Zaragoza, Spam 
^Department of Physics, University of Illinois, Urbana, IL 61801, USA 
5 Laboratory of Elementary Particle Physics, Cornell University, Ithaca, NY 14853, USA 
''School of Physics, University of Edinburgh, King's Buildings, Edinburgh EH9 3JZ, UK 
7 Department of Physics, The Ohio State University, Columbus, OH 43210, USA 

(Dated: October 20, 2010) 

We give the first accurate results for B and B 3 meson masses from lattice QCD including the 
effect of u, d and s sea quarks, and we improve an earlier value for the B c meson mass. By using the 
Highly Improved Staggered Quark action for u/d, s and c quarks and NRQCD for the b quarks, we 
are able to achieve an accuracy in the masses of around 10 MeV. Our results are: tub = 5.291(18) 
GeV, m Bs = 5.363(11) GeV and m Bc = 6.280(10) GeV. Note that all QCD parameters here are 
tuned from other calculations, so these are parameter free tests of QCD against experiment. We also 
give scalar, Bg , and axial vector, B s \, meson masses. We find these to be slightly below threshold 
for decay to BK and B* K respectively. 



I. INTRODUCTION 

B meson physics is one of the critical elements of the 
flavor physics programme. The mass differences between 
the 'heavy' and 'light' eigenstates of the neutral B and 
B s mesons are now known experimentally and can be 
used to precisely constrain the ratio of CKM elements 
|Vtd|/|V( S | if the appropriate theory results have been cal- 
culated with a matching error. A first lattice QCD cal- 
culation of these mixing matrix elements, including the 
effect of u, d and s sea quarks, was given recently and the 
critical quantity £ = fg s \/BbJ JbVBb was obtained to 
3% pQ. Similarly experimental results for the annihila- 
tion of charged B mesons to leptons via a W boson can 
be used to constrain V u b and B semileptonic decays to ir 
or D(*> can be used to constrain V u h and V c b if the ap- 
propriate theory calculations of decay constants or form 
factors are known. Again calculations of these in full lat- 
tice QCD have been done and, for example, the decay 
constant of the B meson is obtained to 7% pQ. 

To improve on these lattice QCD results requires pin- 
ning down and eliminating sources of systematic error 
and testing as stringently as possible that this has been 
done. Here we provide a calculation of the masses of B 
mesons along with their decay constants using b quarks 
in the NRQCD formalism [5] with light quarks in the 
Highly Improved Staggered Quark (HISQ) formalism [3J. 
The HISQ formalism has improved discretisation errors 
compared to the asqtad improved staggered quark for- 
malism used in our previous calculations. Because we 
use NRQCD which can handle hadrons with either sin- 
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gle or multiple b quarks we are able to do an accurate 
calculation of the B meson masses by linking them to 
meson masses in bottomonium. This provides a strong 
test of systematic errors. We are able to handle all of 
the 4 lightest quarks - u/d, s and c - using the HISQ 
formalism and are therefore also able to calculate mass 
differences and decay constant ratios accurately between 
B, B s and B c mesons. 

Section [Tl] outlines how the Lattice QCD calculation 
was done and Section |III| describes the analysis and re- 
sults. Section [TV| compares the results to experiment and 
to other lattice QCD calculations. Section [V] gives our 
conclusion. 



II. LATTICE QCD CALCULATION 

We are concerned here with mesons with one valence 
b quark and a lighter valence anti-quark, either c, s, or I. 
We use the notation I, "light" , to refer to either the u or d 
quark. Everywhere in these calculations mi = m u = m<j, 
but we will correct for the effects of this, along with the 
effects of missing electromagnetic interactions, when we 
compare to experiment. 

The bottom quark moves sufficiently slowly inside 
bound states that it is well described by a non-relativistic 
formulation (NRQCD) on lattices of moderate lattice 
spacing. The lighter partner quark, with balancing mo- 
mentum in a meson with zero total momentum, is moving 
much faster and requires a relativistic formulation. For 
this we use the HISQ formulation, which offers control of 
discretization errors to the level where even the c quark 
can be treated relativistically. 
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Set g a (fm) au 0P m a l sq au 0P m a a sq L/T N conf x N t 

1 6.572 0.1583(13) 0.0097 0.0484 16/48 628 x 2 

2 6.586 0.1595(14) 0.0194 0.0484 16/48 628 x 2 

3 6.760 0.1247(10) 0.005 O05 24/61 507 x 2 

4 6.760 0.1264(11) 0.01 0.05 20/64 589 x 2 

5 7.090 0.0878(7) 0.0062 0.031 28/96 530 x 4 



tion to use for I and s quarks. 

We have shown that the HISQ action can even be used 
for c quarks [3J, but in that case an additional step is 
needed. The 'ordinary' tree-level 0(a 2 ) discretization 
errors coming from the finite difference discretization of 
the covariant derivative are eliminated in both HISQ and 
Asqtad formulations using an additional 3-link 'Naik' 
term. The Naik term corrects errors that would other- 
wise appear at 0(pa) 2 in the quark, and therefore meson, 
dispersion relation of energy vs momentum. A nonper- 
turbative value for the Naik term coefficient (written as 
1 + e) can be derived by studying the dispersion relation 
of the rj c meson and tuning the coefficient until the square 
of the speed of light in this relation is 1. Here we use the 
values of the Naik coefficient determined in this way in 
[10] appropriate to the values of m c a on each ensemble 
|37j . Tuning the Naik coefficient in this way removes all 
discretisation errors from the HISQ action at leading or- 
der in the square of the velocity of the c quark [3J. At 
subleading order in v c there will be 0(a s a 2 ) errors. 

On each configuration of the ensembles in Table [I] we 
have generated and stored random-wall source charm, 
strange and light valence propagators. These were calcu- 
lated in the following way. On the source time-slice to we 
generate a U(l) vector of random numbers r](t ) x >. Then 
we invert to get the HISQ propagator <? HIS Q(x, to): 

g ms ^(x,to) = M-l, V (t ) x , (1) 

For the fine ensemble we use a different random-wall 
at each of four values of to on every configuration. For 
the coarse and very coarse configurations we use two 
random-wall sources per configuration. Although the 
time sources are equally spaced, their position in time 
varies from configuration to configuration through the 
ensemble in a pseudo-random manner to further reduce 
auto-correlations within the ensemble. 

We use the charm, strange and light HISQ quark prop- 
agators from [TO], but note that here we are using a dif- 
ferent definition of the lattice spacing and this affects the 
meson masses in physical units and therefore the tuning 
of the quark masses. We have included some additional 
quark mass values on the very coarse lattices to be able 
to correct for mistuning. We list the HISQ valence pa- 
rameters used in columns 3 to 7 in Table [ill 

The tuning of the HISQ valence c and s masses to 
their correct values on each ensemble is important to 
avoid mistuning effects masquerading as, for example, 
lattice spacing artefacts or sea quark mass dependence. 
We tune the c mass to give the correct r\ c meson mass on 
each ensemble. Here the correct r] c mass has to be ad- 
justed slightly from its experimental value of 2.980 GeV 
to allow for the fact that our lattice QCD calculation 
is happening in a world without electromagnetism and 
without c quarks in the sea and in which we do not allow 
the rj c meson to annihilate to gluons (because we have 
not included such 'disconnected' pieces in our rj c meson 
correlators). These effects are all small and can be esti- 
mated from potential models or perturbation theory. We 



TABLE I: Ensembles (sets) of MILC configurations used with 
gauge coupling /3, size L 3 x T and sea masses (x tadpole pa- 
rameter, uop, taken from the average plaquette) m^ sq and 
m"'. The lattice spacing values in fm are determined us- 
ing the ?7 S meson mass and decay constant [5] and given in 
column 3. Column 7 gives the number of configurations and 
time sources per configuration that we used for calculating 
correlators. On set 5 only half the number were used for light 
quarks. 



A. The gluon configurations 

We make use of the MILC collaborations 's library of 
2 + f -flavor gauge configurations (J. These have two 
degenerate flavors of light sea quarks and one flavor of 
strange sea quark, formulated with the Asqtad action 
[5H?]. The gluon action is Symanzik- improved through 
0(a s a 2 ) except for terms of C(n/a s a 2 ) where n/ is the 
number of sea quarks. In fact these terms [8] are of sim- 
ilar size to the other a s a 2 terms so in practice the gluon 
action has a s a 2 discretisation errors. For this work we 
use five different ensembles at three different lattice spac- 
ings, with a « 0.15,0.11, and 0.09 fm. We refer to these 
as "very coarse", "coarse" and "fine", respectively. The 
configurations have large spatial volumes (> (2.4fm) 3 ). 
Table [I] lists the specific ensembles used in this work. 

In this work we use values of the lattice spacing, a, on 
each ensemble determined using the mass and decay con- 
stant of the pseudoscalar ss meson, the rj s . Although this 
particle is not seen in the real world because of mixing 
with uu and dd which can be prevented on the lattice, 
its properties can be determined from those known from 
experiment of the it and K mesons, as described in [S]. 
Table [I] lists the values obtained in [9] for the ensembles 
we are using here. The values of a are larger on coarse 
lattices than those from the more traditional way of set- 
ting the lattice spacing using the parameter 7*1, but the 
results agree, as they should, in the continuum limit [9j. 

B. HISQ valence quarks 

We use the Highly Improved Staggered Quark (HISQ) 
formulation [3J [TO] for valence charm, strange and light 
quarks. 

The HISQ action further reduces the residual 0(a s a 2 ) 
discretization errors coming from taste-changing effects 
found in the Asqtad formulation. It does this with an 
additional fattening step applied to the gluon field cou- 
pling to the quarks [3J 110) . The errors are reduced by 
about a factor of 3, making HISQ therefore a better ac- 
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Set 






arn c 


1 + e 


am s 


hisq 

am oi 


1 


3.4 


0.8218 


0.85 


0.66 


0.066 


0.0132 




3.6 


0.8218 


0.88 


0.64 


0.08 






3.4(c, / 1) 


0.8218 










2 


3.4 


0.8225 


0.85 


0.66 


0.066 


0.0264 


3 


2.8 


0.8362 


0.65 


0.79 


0.0537 


0.0067 


4 


2.8 


0.8359 


0.66 


0.79 


0.05465 


0.01365 


5 


1.95 


0.8541 


0.43 


0.885 


0.0366 


0.00705 



TABLE II: Parameters for the valence quarks, amt is the b 
quark mass in NRQCD, and uol is the tadpole-improvement 
parameter used there |llj . We use stability parameter [11] 
n = 4 in NRQCD everywhere. Columns 4, 6 and 7 give the 
charm, strange and light bare quark masses for the HISQ 
action. 1 + e is the coefficient of the Naik term in the charm 
case [3]. On set 1 we give additional values of mt and m 3 that 
were used for tuning purposes as described in the text. We 
also used alternative d coefficients for the arrib = 3.4 case, 
specifically ci = eg = 1.36 and C5 = 1.21, again as described 
in the text. 



take the physical value of m^ c appropriate to our cal- 
culations to be 2.985(3) GeV [IS], incorporating a shift 
of 2.6 MeV for electromagnetic effects and 2.4 MeV for 
annihilation effects. The effect of c quarks in the sea is 
negligible; this will be discussed in the next section. 

We tune the s quark mass from the value of the mass 
of the rj s meson. This is a fictitious pseudoscalar meson 
made of s quarks, which is not allowed in our lattice cal- 
culation to mix, by annihilation, with uu and dd mesons. 
We cannot therefore fix its mass from experiment but 
must do within a lattice QCD calculation, in which we 
extrapolate results for 7T, K and T) s meson masses and 
decay constants simultaneously to the continuum limit 
and physical point for the I and s quarks. This was done 
in [S] and the value M Va = 0.6858(40) GeV was obtained, 
as the value appropriate to this lattice QCD world. 

The light quark valence masses are taken to match ap- 
proximately the light quark masses used in the sea. The 
way that this is done is to take the ratio of the light va- 
lence HISQ mass to the appropriate HISQ strange mass 
value to be the same as the ratio of the sea light quark 
mass to its appropriate s mass value (the Ugp factors can- 
cel in that ratio). This can only be done approximately 
because the sea 'strange' quark masses are not very well 
tuned in some cases [4] and the correct value for s is not 
known very accurately. However, sea quark mass depen- 
dence is only a very small effect for everything calculated 
here, so this is not a big issue. 



C. NRQCD b quarks 

One can estimate the velocity of the quarks in bot- 
tomonium mesons by comparing radial excitation ener- 
gies to masses. This shows that in the T the b is very 
non-relativistic with v\ ~ 0.1 (in units of c 2 ). By com- 
parison, in charmonium v 2 c sa 0.3. In a mixed system 



with a lighter quark the b quark is even slower than in 
bottomonium. Consider that in a B c meson the reduced 
quark mass is roughly 1.5 times that of the bb system 
and 4.5 times that of a cc system. For a constant mean 
kinetic energy |13j across all three systems, we then ex- 
pect v\ 1=3 0.04 and v\ s» 0.35 in the B c . For a b quark 
combined with an even lighter quark, v\ will get even 
smaller. Assuming a mean momentum of Aqcd for the 
b quark inside a B or B s would give v\ ~ 0.01. An 
NRQCD approach is then well justified for the b quark 
in B c , B s and B\ mesons when combined with a relativis- 
tic approach for the lighter quarks. 

As we have used a random-wall source for the HISQ 
propagators, it is critical that we initialize the NRQCD b 
propagators b with the same random-wall function r)(t ) x i 
as we used for the HISQ propagators. This is non-trivial 
because the HISQ propagators have one only Dirac com- 
ponent, whereas NRQCD propagators have explicit spin. 
The spin information for the HISQ propagators is tied up 
with the positions of source and sink, however the source 
site information is lost for a random wall source once the 
propagator has been calculated. If we had access to a 
HISQ propagator from all source points to all sink points 
we could reconstruct the full 4x4 spin structure by multi- 
plying by Q(x)Q^ (y) where Q is the staggering operator, 
given as a product of Dirac gamma matrices as: 

tt(*) = 7r 7 2 X2 7 3 X3 7 t Xt - ( 2 ) 

Here we cannot apply the source fi to the HISQ propa- 
gators, but instead we can apply it as a source for the 
NRQCD propagators. This then effectively 'undoes' the 
staggering transformation and gives a naive quark source 
that can be combined with a b quark source, but it is done 
after the staggered propagators have been made. The b- 
quark source then has to have 4 spin components, rather 
than the usual two for NRQCD because we cannot sepa- 
rate the upper and lower components of the naive quark 
source after the fact. Our b quark source is constructed 
from Q(x) multiplying the appropriate random noise at 
each site, r](to) x >, updating the standard method of com- 
bining staggered quarks with other formalisms that have 
explicit spin components derived in |14j . 

A further issue is that we must project onto pseu- 
doscalar or vector heavy-light mesons directly at the 
source of the b quark propagator rather than combin- 
ing appropriate spin components at the end. We then 
have to calculate separate b quark propagators for each 
of the pseudoscalar and vector mesons. 

Finally, to enhance our ability to isolate the ground- 
state behavior, we smear the b propagator source with a 
smearing function, S, of various functional forms and dif- 
fering radii, r,. Combining all of these factors we there- 
fore, on timeslice to, initialize the NRQCD propagator 
as: 

G NRQCD (X)fo) = £ s(lx _ ^ ;r . )7?x ,( to)r0(a ,' )) (3) 
x' 

where T is an element of the Dirac algebra chosen to 
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project out either a pseudoscalar or a vector heavy-light 
meson state. 

On subsequent timeslices we evolve the NRQCD prop- 
agator recursively in the standard way [TS]. Note that 
here, because upper and lower components do not mix 
in NRQCD, the upper and lower halves of our b-quark 
source are simply evolved separately with the same evolu- 
tion equation. Our b-quark propagators then have only 2 
spin components at the sink end. The evolution equation 



Gi(x,t+1) = 
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5H = - 
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_B and i? are improved versions of the naive lattice 
chromo-electric and chromo-magnetic fields, E and B. 
All the gauge fields appearing are tadpole-improved by 
dividing by a tadpole factor, uql, derived from the mean 
link in Landau gauge. 

The equations given above represent the standard 
NRQCD action, used in many previous lattice QCD cal- 
culations (for example [15] and [H]), correct through v^. 
The largest source of remaining systematic error from 
this action is from radiative corrections to the coefficients 
of the v$ terms required to match full QCD through 
0{a s vf) (the v 2 term is tuned nonperturbatively in fixing 
the b quark mass as described below) . Here we generally 
use the tree-level values of c, = 1 for the constants as 
we have done before. However, we have also done some 
calculations on set 1 using values of C\ , C5 and Cq that in- 
clude radiative corrections, to be able to gauge the size of 
the systematic error from these terms. How the radiative 
corrections are calculated and further tests of them will 
be described elsewhere [TBI IT7] . The values we use on set 

1 are those appropriate to a value of m^a of 3.4 and a s in 
the V scheme appropriate to the lattice spacing of that 
ensemble. These are c\ — c§ = 1.36 and C5 = 1.21 [16] . 

We list the NRQCD valence b masses used in column 

2 of Table |TT] along with the uol parameters. We have 
used two different masses on set 1, again to test for sys- 
tematics from quark mass tuning. Since NRQCD quarks 
propagate in one direction in time only we improve statis- 
tics by generating propagators both forwards in time (for 



T/2 time units) and backwards in time from each light 
quark source. 

The b quark mass is tuned by determining T and rjj, 
meson masses [T7]. Because the zero of energy has been 
shifted in NRQCD we cannot determine meson masses 
directly from their energy at zero momentum as is done 
in a relativistic formulation. Instead we must calculate 
the 'kinetic mass' of a meson, M, which appears in the 
relationship between E(p) and p 2 : 



E{p) = E{0) + vV + M 2 - M. 



(7) 



We are able to do this very precisely using random wall 
sources [5] ■ We determine the kinetic mass of both the T 
and the rjb mesons and use the spin-average of the two, 
i.e. 



M b - b 



3M T + M v 



(8) 



to tune the b quark mass. The reason that we use the 
spin-average is to avoid systematic errors from the terms 
in the NRQCD action that give rise to spin-splittings. 
These terms are only included at leading (vf) order in 
the action above (equation [6]), and we know that there 
are sizeable discretisation errors in the hyperfine splitting 
between the T and f]t as a result [TT] . These errors also 
have a small effect on the kinetic masses [T7] and we take 
the spin-average to remove them. 

The experimental value for M b -t that we tune to must 
be adjusted, as described above in the r\ c case, to be 
the value appropriate to our lattice QCD world which is 
missing some elements of the real world. We estimate 
the effect of the electromagnetic attraction of the b and 
b to be 1.6 MeV from a potential model. The absence 
of electromagnetism shifts both My and M Vb upwards. 
The shift from not allowing the r/f, to annihilate to gluons 
we take to be the same as for the r\ c at 2.4 MeV. In 
addition we must consider the absence of c quarks in the 
sea. This has only a very small effect on charmonium 
states [12] but a somewhat larger effect on bottomonium, 
because these states are more sensitive to exchange of 
higher momentum gluons that can generate a cc pair. 
The effect can be estimated perturbatively from that of 
a massive quark loop in the gluon propagator that gives 
rise to the heavy quark potential [T2]. It gives rise to 
a correction to the potential which is proportional to a 
delta function at the origin: 



V(r) = - 



C 



-c f 



a 
10m 



(9) 



This correction is very similar to the hyperfine potential 
in bottomonium, except that it is not spin-dependent 
and the mass that appears in the denominator is m c and 
not mb- The hyperfine potential induces the mass differ- 
ence between My and M nb of 69 MeV [18] . The potential 
above has a coefficient 14 times smaller, where this factor 
is given by [807r/(3a s )] x (m c /m(,) 2 and we take mb/m c 
= 4.51 from [19] ■ We therefore expect a shift of around 5 
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MeV to both T and 77^. Again the absence of c quarks in 
the sea shifts the mass values upwards. The experimen- 
tal % mass is now 9.391(3) GeV [5D] and the T, 9.460 
GeV, giving a spin-average of 9.443(1) GeV. Applying 
the shifts above with 50% errors, we find the appropriate 
value for M b r for us to tune to is 9.450(4) (1) GeV, where 
the first error is from applying the shifts and the second 
is the experimental error. 

D. HISQ-NRQCD two-point functions 

At the sink end of the propagators we must tie a two- 
component NRQCD propagator with the one-component 
HISQ propagator. We first convert the HISQ propagator 
back into a multi-spin object by multiplying by Cl(x) at 
each sink site: 

G HIS( 3(x, t) ab = g HISQ (x, i )n(x, t) ab , (10) 

remembering that the W factor that would normally be 
at the source end has now been included in the NRQCD 
propagators. 

Finally we generate the B meson correlator by com- 
bining the HISQ propagator and NRQCD propagators at 
the sink timeslice t, again with smearing functions and 
appropriate Dirac structure: 

Cr(t - hh = E G HISQt (x, t)TS(\x - x'\ ; r 3 ) (11) 

xG NRQCD (Xjt) 

Note that this is a sum over 4 spin components at the 
source and 2 at the sink, so T is either a 2 x 2 unit matrix 
for the pseudoscalar meson or a Pauli spin matrix for a 
vector meson. 



E. Extracting physical masses 
We do simultaneous constrained fits [3T] to the form 

«cxp 

C meson (i,j,t ;t) = Y. h *- kh lk e ~ Ek{t ~ ta) ( 12 ) 
fe=i 

k' = l 

where i and j index the 3x3 matrix of source and sink 
smearing functions. The second term fits the oscillating 
"parity-partner" states that appear in most staggered 
meson correlators. Our correlation functions cover the 
range of t — t Q values from zero to T/2, although we do 
not fit all the way to t — to = 0. Instead we start at 
tmin = 2 — 4 for B s and Bi fits and 6 — 8 for B c fits to 
reduce the effect of excited states. 

We constrain the parameters of the fit with prior val- 
ues and widths, which are fed into the augmented x 2 



function that the fit minimizes. These priors represent 
very general information about mass splittings and am- 
plitudes. The prior value for the ground state mass is 
simply taken from an effective mass plot, with the prior 
width taken to be a factor of 1.5 from this value. The 
mass parameters for the higher mass states enter the fit 
as the logarithm of the mass difference with the state 
immediately below, so that the mass difference is posi- 
tive and higher mass states remain, by definition, higher 
mass. We take a prior value on these mass differences, 
both for normal states and oscillating states, to be ps 600 
MeV, converted to lattice units for each fit. The prior 
width on the mass differences is taken to be a factor of 
2. The prior on the mass of the lowest oscillating state 
is taken to be ~ 400 MeV above the ground state, with 
a prior width of a factor of 1.5. We use 0.1 ±1.0 for 
the prior value and width for all amplitudes. Here the 
prior value of 0.1 is simply to provide a non-zero start- 
ing point for the fit. The width of 1.0 can be uniform 
across correlators of different smearing functions because 
we normalise them all so that S(x) 2 — 1 across a 
timeslice. 

In this way, we are able to obtain high-confidence fits 
which are stable, both in the central value for the ground 
state mass and amplitude and their errors, with respect 
to varying the number of exponential functions included, 
N cxp . We take our results from fits with N exp = 5, since 
all our fits are stable by this point. 

Where possible (sets 2-4) we simultaneously fit mul- 
tiple light valence channels. That is, by fitting simulta- 
neously Cb b and Cb u we can eliminate correlated errors 
from estimations of mass differences such as Mb, — ■ 

One important issue with B meson correlators is their 
exponentially falling signal/noise ratio, which means that 
the statistical accuracy that can be obtained on masses 
and decay constants is not as high as that of lighter 
mesons, for example D/D s [12]. The variance of the B s 
correlator, for example, contains bbss propagators and 
can rearrange them into an r)b and an rj s . Thus the noise 
(square root of the variance) falls exponentially with a 
lowest energy (E nb +m Vs )/2 at large times while the sig- 
nal falls with a lowest energy Eg s . This means that the 
signal/noise degrades exponentially with a physical en- 
ergy which is the mass difference between the B s and 
(M m +M na )/2 (330 MeV). This is illustrated in Figure [I] 
where we explicitly compare the effective mass of the Bj 
correlator and the effective mass of its statistical error, 
and show that the 'mass in the noise' is as expected. This 
physical mass difference cannot be altered but if we use 
smearing functions, as we have done here, it is possible 
to extract the ground state B s mass from early t — to 
values, where the noise is less of an issue. 

As discussed earlier, the zero of energy is changed in 
the NRQCD formulation so the energy parameters E\. 
and E'ui include an energy shift for which we must correct 
before comparing to experimental values. We can do this 
by comparing the B meson state of interest (containing 1 
b quark) to a reference state, which can also be calculated 
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FIG. 1: The effective mass of the signal (left) and noise (right) for B s correlators on the fine ensemble (set 5). The open squares 
are from correlators with local sources and sinks and the bursts from correlators with smeared sources and sinks. The effective 
mass is obtained from the natural logarithm of the ratio of the correlator (or its error) on successive time slices. The blue line 
on the left plot corresponds to the fitted energy of the B s meson and the red line on the right plot to one half the sum of that 
for the rjb and rj 3 mesons (at the same ami and am 3 ). This figure should be compared to Figure 3 in [12] . 



with NRQCD b quarks on the lattice and whose mass is 
known measured experimentally. That is, 

M B = (E B - -E^a- 1 + -M ref , (13) 
n n 

where E TC f is calculated on the same lattice ensemble, 
and -Mr e f comes from experiment (adjusted if necessary 
for the absence of electromagnetism etc. from our cal- 
culation), n is the number of b quarks in the reference 
state. The reference state can also be a linear combina- 
tion of states, such as the spin-average of bottomonium 
states that we will use below. 

To minimize the contribution of the 0.8% uncertainty 
on a -1 to the overall uncertainty in Mb, it is important 
to choose a reference state that makes the quantity (E B — 
E xe f/ii) as small as possible. We will sometimes do this 
below by subtracting the masses of additional reference 
states, for example ones made purely of c quarks for the 
B r case. 



III. RESULTS ON THE SPECTRUM 

A. B s meson mass 

To determine Mb s we follow the strategy described 
earlier, using the spin-average of bb states as a reference, 
and calculate: 



Ar =\E r - iRx ) a- 1 



latt 



From this we can reconstruct Mb s using 

1 

2 J "b&,phys 



M Bs ,Utt = A Bs + rMit . 



(14) 



(15) 



on the same configurations as used for calculating the B s 
meson energies. M bb is used to tune the b quark mass, 
as discussed earlier, and M bbphys is the physical value 
taken from experiment, but adjusted (to 9.450 GeV) for 
the lattice QCD world (missing electromagnetism, r^, an- 
nihilation and charm-in-the-sea). To compare our results 
for Mg st i a tt to experiment we have to add corrections to 
put back in missing electromagnetism and charm-in-the- 
sea effects. These corrections are negligible, however, as 
we will discuss below. 

Table Mil lists all our fitted values needed for determi- 
nation of the B s meson mass. Note that the error on 
the fitted B s meson energy is larger than any of the er- 
rors on the other fitted energies. This is because of the 
signal/noise problem in the B s correlator discussed ear- 
lier. The other correlators used here do not have that 
problem and the fits give much more precise results for 
ground state masses. Details of these other fits are given 
elsewhere [9]. 

In Figures [2] and [3] we show how A b s varies with the 
square of the rj s mass and the spin-averaged bb mass from 
our results on set 1. These results allow us to correct for, 
and estimate the errors from, mistiming b and s masses. 
The lines are simple linear fits in M bb and M? . The slope 
of A B „ against M 2 is 0.19, in good agreement with that 
expected from the experimental data comparing B s and 
B. The slope against M bb is very small because the b 
quark mass effects naively cancel in Ab 3 ■ However some 
residual dependence remains and gives a slope of 0.017, 
somewhat smaller than the experimental result of 0.033 
obtained over a much larger mass range from comparing 
B s and D s . 



Table IV gives the values of Ab 3 , adjusted for mistun- 



Here E bb is the spin-average of T and % energies at zero 
momentum calculated with the same NRQCD action and 



ing by using the slopes given above and the mismatch of 
M bb and on each ensemble compared to the phys- 
ical values (9 S .450 GeV and (0.6858 GeV) 2 respectively). 
We take a 50% error on any shift applied for mistuning. 
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Set 


arrib 




(lEi T 


0,7X1 s 


aM r]s 




aAs VP 


aA° + -°" 


aAj + - 1_ 


1 


3.4 


7.260(9) 


0.27843(8) 


0.066 


0.52524(36) 


0.6409(11) 


0.0343(11) 


0.310(11) 


0.308(14) 




3.4 






0.080 


0.57828(34) 


0.6539(10) 


0.0349(9) 


0.317(11) 


0.309(12) 




3.6 


7.688(5) 


0.27662(7) 


0.066 




0.6466(13) 


0.0324(14) 


0.300(15) 


0.307(16) 




3.6 






0.080 




0.6604(9) 


0.0332(9) 


0.315(11) 


0.308(11) 




3.4( Cl £ 1) 


7.248(4) 


0.28048(7) 














2 


3.4 


7.261(9) 


0.27902(7) 


0.066 


0.52458(35) 


0.6417(14) 


0.0359(21) 


0.299(17) 


0.316(19) 


3 


2.8 


5.996(8) 


0.28538(3) 


0.0537 


0.43118(18) 


0.5470(15) 


0.0287(19) 


0.215(17) 


0.254(6) 


4 


2.8 


5.992(5) 


0.28465(6) 


0.05465 


0.43675(24) 


0.5527(16) 


0.0261(27) 


0.253(8) 


0.249(16) 


5 


1.95 


4.288(10) 


0.25985(5) 


0.0366 


0.30675(12) 


0.4172(10) 


0.0189(12) 


0.1708(48) 


0.166(11) 



TABLE III: Results for energies and kinetic masses needed for the determination of the mass of the B s meson. All the results 
are in units of the lattice spacing. For each set we list the valence b quark mass and its associated kinetic mass and energy for 
the spin-average of T and rjb states. We also give the valence s quark mass and its associated rj s meson mass. These values 
are also given in [9]. Where we have used multiple 6 and s masses on set 1, we give the r] 3 and E bb values only once to avoid 
confusion. In column 7, we give the fitted energy of the B a meson (i.e. Eq from fits to th e form given in equation [l2|. In 



hid 



column 8 we give the hyperfine splitting, A^ vp — E(B*) — E(B a ), discussed in section 
but includes some additional values on set 1 that we use for studying systematic errors 
line 1 is different from that in 22 although consistent with it 

alone, rather than from a full fit including Bi and B* , to be in keeping with the other B 



This column is largely from [22] 
Note also that the value for set 1 on 
Here we use a value from a fit to the B 3 and B* correlators 

values given on set 1. Columns 9 



and 10 give the values of mass differences between scalar and pseudoscalar and between axial vector and vector respectively, 
discussed in section UlI El 
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FIG. 2: Results from set 1 for the mass of the B s meson 
(specifically the difference between that mass and one half of 
the spin-averaged mass of T and rjb) as a function of the square 
of the r/s meson mass, acting as a proxy for the strange quark 
mass. The errors are statistical only, since lattice spacing 
errors affect all the points together. The lines are fits to the 
results allowing linear terms in M% and M bb . Here the lines 
join points for a fixed b quark mass. See Figure [3] for the 
equivalent as a function of M bb . 



Errors from mistuning are smaller than the statistical er- 
rors except on sets 1 and 2. Table [TV] also gives the error 
from the uncertainty in the lattice spacing from Table [T] 
The error is smaller by a factor of 2 than the naive result 
of multiplying A^, by the percentage error in a . The 
reason is that changing the lattice spacing requires the 
masses to be retuned and this affects A in the opposite 



FIG. 3: Results from set 1 for the mass of the B s meson 
(specifically the difference between that mass and one half of 
the spin-averaged mass of T and r/b) as a function of the spin- 
averaged kinetic mass of the T and rjb, acting as a proxy for 
the b quark mass. The errors are statistical only, since lattice 
spacing errors affect all the points together. The lines are fits 
to the results allowing linear terms in M? g and M bb . Here the 
lines join points for a fixed s quark mass. See Figure [2] for the 
equivalent as a function of M^ g . 

direction. 

The resulting errors on the tuned values for Ab s are 
typically less than 1%, around 4 MeV. Within these un- 
certainties we are not able to distinguish any dependence 
on sea quark masses or the lattice spacing. Sea quark 
mass effects are expected to be very small, because the B s 
is a gold-plated particle and has no valence light quarks. 
The lattice spacing dependence depends on the quantity 
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Set 


A Ss (GeV) 5 Xl 5x s 


1 

2 


0.6392(16)(28)(25) 0.15 -0.08 
0.6382(17)(85)(29) 0.34 -0.08 


3 
4 


0.6401(24)(02)(26) 0.09 0.29 
0.6433(25)(13)(29) 0.22 0.29 


5 


0.6417(22)(20)(26) 0.18 0.09 



TABLE IV: Results for A_b s (the mass difference between 
the B a meson and the spin average of T and r/i, masses) on 
different ensembles after tuning to the correct valence b and s 
quark masses. The 3 errors listed are statistical, tuning and 
from the uncertainty in the lattice spacing. Columns 3 and 
4 give Sxi and Sx s , the fractional mistuning of the sea quark 
masses in units of the s quark mass, as defined in the text. 



the coefficients of discretisation corrections can be arrib 
dependent in NRQCD and therefore our discretisation er- 
rors can be amb dependent. We need to allow for a mild 
nonsingular dependence (i.e. appropriate to the values of 
amb that we are using) in our fits, so that the systematic 
error from this can appear in our final results. Since 
any smooth function can be expanded over a limited 
range using a polynomial, we simply allow for linear and 
quadratic terms in the variable Sx m = {ra^a — 2.7)/1.45. 
The factors 2.7 and 1.45 are chosen so that 8x m changes 
from -0.5 on the fine lattices to +0.5 on the very coarse 
lattices as m^a covers the range that we have used. 

We therefore fit our results for to the following 
form: 



chosen to fix the lattice spacing. Earlier reporting of 
these results [53], using the variable r\ to set the scale, 
did show visible lattice spacing dependence. Here it ap- 
pears, perhaps not surprisingly, as if Ab s has the same 
discretisation errors as the i] s used to the fix the scale. 

The tuned results from Table IIVI are used to recon- 
struct Ms s ,iatt (using equation 15) and this is plotted 



A Bg (a,Sxi,Sx s ) = A Bs! p hys [l 



(16) 



in Figure [4] against the square ofthe lattice spacing. In 
order to quote a physical value that can be compared 
to experiment we need to fit our results as a function of 
lattice spacing and sea quark mass so that systematic er- 
rors from such dependence can be fed into our final error. 
The sea quark mass dependence we take to be a simple 
polynomial form in the variables 8x s and Sxi , defined by 
Sx q = (m giSea - m g , S ea,phys)/"Vsea,phys- These variables 
were used in [TJ] but must be adjusted here consistently 
for the change in definition of the lattice spacing and the 
new values are given in Table |IV| Any sea quark mass 
dependence identified in our fit can be extrapolated to 
the physical point where Sxi = dx s = 0, and our errors 
allow for dependence not resolved by our fit. 

The lattice spacing dependence is a trickier issue in 
NRQCD because we cannot extrapolate naively to a = 0. 
What we need to do is to fit the lattice spacing depen- 
dence and assess, using information from the fit, how 
much of the dependence is physical and how much un- 
physical, and allow for both in the final error. Physical 
dependence on the lattice spacing will arise from discreti- 
sation errors in the gluon and sea quark actions, and in 
the light valence quark (HISQ) action. We expect this 
dependence to be 0(a s a 2 ) at leading order as discussed 
earlier. 

The NRQCD action also has discretisation errors. 
These are corrected at tree level by the terms with co- 
efficients C5 and cq in equation [6j Beyond tree level C5 
and eg have an expansion in powers of a s , required for 
NRQCD to match QCD at that order, whose coefficients 
depend on ami,. This dependence will typically be mild 
for large ami but become singular as ami, — > 0. This 
has been explicitly checked for the C5 and eg coefficients 
for a slightly different action in [24] and results have also 
been derived for this action [If}] [T7] . They show almost 
no ami, dependence for ami, > 1. In general, however, 



+ ^ Cj(Aa) 2j (l + Cj b Sx m + C-jbbiSXmf) 

+ 2biSx t (l + Cl (Aa) 2 ) 

+ 2b s Sx s {l + c s (Aa) 2 ) 

+ Ab u (8xi) 2 + 2b ls Sx l 6x s + b ss (5x s ) 2 }. 

We take the prior on AB 3 ,phys to be 0.6(2). The priors 
on the sea quark mass dependence, bi and b s , are taken 
to be 0.00(7). Sea quark mass effects are suppressed 
by a factor of 3 over valence mass effects and here va- 
lence effects correspond to a slope in quark mass of less 
than 0.2. The priors for the b parameters correspond- 
ing to the quadratic sea mass dependence are then set to 
(0.2) 2 /3 i.e. 0.000(13). We take the scale ofthe physical 
a-dependence to be the scale of A = 400 MeV, since we 
expect it to be set by typical internal meson momenta in 
QCD. The coefficients of the quadratic a dependence, c\, 
a and c s , should be 0(a s ) so we take priors of 0.0(5). 
For C2, and the amb dependence of the discretisation er- 
rors, Cjb and Cjbb, we take a very conservative prior of 
0(1). 

The result from the fit is A BsiPhys = 0.638(6) GeV. 
The fit sees no dependence on lattice spacing, amb or sea 
quark masses but our final error allows for all of these. 
The fit, and its error, is robust under changes in the 
number of fit parameters, for exam ple , including or not 



including the a terms in equation 16 It is also robust 



under changes the prior widths. For example, doubling 
the prior width on the lattice spacing or sea quark mass 
dependence changes the final result by less than 1 MeV. 

In the error budget in Table [V] we separate the 6 MeV 
error into component parts coming from the errors on the 
original data points (statistics, tuning and uncertainty 
in the lattice spacing) and the errors coming from the 
lattice spacing and sea quark mass dependence of the 
results, using the method described in |12) . The error on 
the original data points dominates. 

To reconstruct Mb 3 ,\&tt we must add 9.450/2 GeV to 



As^phys as in equation 15 This gives Mb s in the lattice 
world with no electromagnetism or c quarks in the sea. 
The latter effect should be negligible for the B s since it 
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FIG. 4: Results for the mass of the B s meson tuned to the cor- 
rect a and 6 quark masses on each ensemble, plotted against 
the square of the lattice spacing. The errors on the points in- 
clude statistical and tuning errors and the uncertainty in the 
lattice spacing. The dark shaded band is our physical result, 
allowing for sea quark mass and lattice spacing dependence as 
described in the text. The width of the light shaded band re- 
flects our full error as given in the error budget, Table M The 
black star is the experimental result [20], offset from a — 
for clarity. 



is a much larger particle than the T or r\b and therefore 
much less sensitive to the gluon exchange that could cre- 
ate a cc pair. The effect of electromagnetism can be esti- 
mated following |12] . There we gave a phenomenological 
formula for electromagnetic and m u /md mass difference 
effects in heavy-light mesons: 

M(Q, q) = M sim (Q,q)+Ae q e Q +Be 2 q +C(m q -mi). (17) 

M s ; m is the mass of the meson in the absence of elec- 
tromagnetism and with m u = m,d- From experimental 
charged and neutral B and D meson masses we deter- 
mined A sa 4 MeV and B w 3 MeV. For the B s then this 
formula gives a shift between M(Q, q) and M sim (Q, q) as 
a result of electromagnetism of -0.1 MeV, a very small 
effect. We make no correction for this, but add an error 
for it to our error budget. 

Additional systematic errors that must be added in to 
the error budget are those from relativistic corrections 
that are not included in our NRQCD action. These er- 
rors affect results at all lattice spacings equally and so 
cannot be estimated from our results as we have done for 
the discretisation errors. We must consider the effect of 
relativistic corrections on both the B s mass and on the T 
and r\b masses because they both appear in Ab s . In fact 
we expect relativistic corrections to have a bigger effect in 
bottomonium than on the B s . Our NRQCD action is cor- 
rect through 0(vf) for bottomonium and so the largest 
missing terms are at a s v^ and vf. We expect the typical 
energy shift of a spin-independent v% term to be ~ 50 
MeV (0.1 x 500 MeV), so a s vf corrections could give rise 



Error 


M Bs M Bo ,hh 


M Bo ,hs 


Stats/tuning/uncty in a 


5.5 


2.9 


14.0 


Lattice spacing dependence 


0.5 


2.9 


8.0 


Tig, sea dependence 


3.0 


1.0 


4.0 


spin-ind. NRQCD systs. 


8.0 


6.0 


7.5 


spin-dep. NRQCD systs. 


3.5 


4.0 


1.0 


uncty in M Vs 


1.0 


- 


2.3 


em, annihiln, c sea in bb 


2.5 


2.5* 


0.0 


em, annihiln, c sea in cc 




1.5* 


0.2 


em effects in B s or B c 


0.1 


1.0* 


1.0 


em effects in D a 






0.7 


finite volume 


0.0 


0.0 


0.0 


Total (MeV) 


11 


9.5 


19 



TABLE V: Full error budget for B s and B c meson masses in 
MeV. See the text for a fuller description of each error. The 
total errors are obtained by adding the individual errors in 
quadrature except in the case of the starred errors. They are 
correlated and this must be taken into account as described 
in the text before being squared and accumulated into the 
square of the total. 



to 15 MeV shifts in M b T. Similarly a spin- independent 
vf correction could give rise to an energy shift of ~ 5 
MeV. Adding these two in quadrature and dividing by 
2 gives an estimate of the systematic error in Mb s from 
relativistic corrections to M bb of 8 MeV. For B s , the ap- 
propriate power counting for relativistic corrections is in 
Vb ~ A/m;,. Our NRQCD action already includes high 
order terms in A/ nib at tree level and so there are no 
significant tree level errors for the B s . The leading error 
is at ct a Vb in missing radiative corrections to C4, the co- 
efficient that multiplies the a ■ B term giving rise to the 
hyperfine splitting. As we will discuss in section III D 



we have plenty of evidence that errors coming from this 
term are small, at most 10% of the hyperfine splitting 
itself. Since this error would vanish for the spin-average 
of the masses of the B s and the B*, which is 3/4 of the 
hyperfine splitting above the B s mass, we take the error 
in the B s mass to be 3/4 of the error in the hyperfine 
splitting, 3.5 MeV. 



Errors from the uncertainties in M„ 



and M,t that we 



use to tune the s and b quark masses can be estimated 
from the slopes in Figures [2] and [3j The 4 MeV uncer- 
tainty in M Vs feeds into a 1 MeV uncertainty in A b s and 
therefore Mb, ■ The uncertainty in A b s from the 5 MeV 
uncertainty (simply adding the statistical and systematic 
errors) in M bb is very small - less than 0.1 MeV because 
of the cancellation of the b quark mass inside A b s ■ How- 
ever, the uncertainty reappears when we reconstruct Mb s 
by adding M bb /2 to A# 3 . This then gives a sizeable 2.5 
MeV error in Mb s ■ 

Finite volume errors are expected from chiral pertur- 
bation theory to be negligible for the masses of mesons 
containing heavy quarks on volumes exceeding (2.4fm) 3 , 
that we are using here. 

The full error budget is given in Table |Vj The system- 
atic errors discussed above, added in quadrature, give 9 
MeV, dominating the 6 MeV errors coming from the sta- 
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tistical errors of the data and its lattice spacing and sea 
quark mass dependence. The final result is then M Bs = 
5.363(6) (9) GeV. Figure [4] shows a dark shaded band for 
the first error and a lighter shaded band to encompass the 
full error, adding 6 MeV and 9 MeV in quadrature to give 
11 MeV. To reduce the full error will require improve- 
ments to the NRQCD action, currently underway. The 
experimental result for the B s mass is 5.3663(6) GeV. 



B. B c mass 

For the B c meson mass we could use exactly the same 
procedure as for the B s . However, there is a better 
method, in which we subtract in addition the mass of 
a charmonium reference state, the rj Cl to reduce the en- 
ergy difference calculated on the lattice to a very small 
value: 
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9.2 9.4 
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9/1 
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E Be -\{E bl + M Vc 



latt 



(18) 



FIG. 5: Results from set 1 for the mass of the B c meson 
(specifically the difference between that mass and one half of 
the spin-averaged mass of T and rjb added to the mass of the 
r/ c meson) as a function of the spin-averaged kinetic mass of 
the T and r/t, acting as a proxy for the b quark mass. The 
errors are statistical only, since lattice spacing errors affect all 
the points together. The lines are fits to the results allowing 



We call this the "heavy-heavy" (hh) subtraction method. linear terms in M ^ and M bb- Here the lines join points for 
Here M„ is the value of the r] c mass calculated on the a fixed c 1 uark mass ' See Fi S ure © for the equivalent as a 



lattice and M„ , p hys 1S its value from experiment appro- 
priately adjusted for the lattice QCD world, as described 
earlier. For the HISQ quarks that we use for c the en- 
ergy obtained from fits to charmonium correlators at zero 
momentum is the charmonium mass, so there is no issue 
with the zero of energy. We simply use the additional 
charmonium subtraction here to reduce errors from the 
uncertainty in the lattice spacing. 

We will compare results of this to a second method: 



function of M r)c . The third point at lower left gives results 
for a ^ 1; it is not included in the fit. 



M Ba = (E Bc - (E Bs +M Ds )) 



latt 



(M 



B s ,phys 



M 



D s ,physJ ? 



(19) 



that we call the "heavy-strange" (hs) subtraction 
method. Here we are using the B s meson to cancel the 
NRQCD shift of the zero of energy in the B c . Again the 
subtraction of the D s meson mass, calculated with HISQ 
c and s quarks is simply to reduce lattice spacing errors 
from the mass difference. 

We first discuss results from the hh method. The B r 
energies and M Vc masses are given in Table 



0.104 
0.102 - 
0.1 

S 0.098 - 
0.096 
0.094 
0.092 - 
0.09 



2.82 2.84 



2.9 2.92 2.94 



VI 



and the bb 

energies and masses, already used in the determination 



of the B s mass, are given in Table III As before, we have 
to tune quark masses on each ensemble to their correct 
value. We show in Figure [5] how the splitting 



A 



B c .hh 



Er 



(E b - b + M v J 



(20) 



latt 



depends on M bb - from our results on set 1 at two values 
of am b and two values of am c . We see that the slope is 
very small, 0.014, because very little b quark mass de- 
pendence is left after cancellation in this mass difference. 



FIG. 6: Results from set 1 for the mass of the B c meson 
(specifically the difference between that mass and one half of 
the spin-averaged mass of T and ryt, added to the mass of the 
r/ c meson) as a function of the mass of the r\ c meson, acting as 
a proxy for the c quark mass. The errors are statistical only, 
since lattice spacing errors affect all the points together. The 
lines are fits to the results allowing linear terms in M Vc and 
A/jj. Here the lines join points for a fixed 6 quark mass. See 
Figure [5] for the equivalent as a function of M^. The third 
point at lower left gives results for a ^= 1; it is not included 
in the fit. 
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An estimate can be derived for the expected slope by 
comparing results for the b quark mass set to the value 
of the c quark mass (when Ag^/j becomes -3/8 times 
the charmonium hyperfine splitting). This gives a slope 
of 0.016 over a much wider range. Figure [5] also shows 
the value of the mass difference for the case where we use 
an NRQCD action with c\, c 5 and Cq set to the values 
including 0(a s ) radiative corrections appropriate for set 
1. We see that this makes negligible difference. 

Figure [6] shows the results as a function of M Vc . The 
slope here is very small but in the opposite direction to 
that for the dependence on M bb -. The value of the slope is 
-0.004. Based on the arguments above we would expect a 
slope of opposite sign but about 60% that of the 6-quark 
mass dependence. However, as stated above, this esti- 
mate is made over a much larger range than the Figure. 



We can use the slope against M bb and against M„ c to 
correct for the slight mistunings of the b quark and the c 
quark that we have on some ensembles. Even though the 
shifts from M bb and M ric dependence are very small they 
are not negligible. This is because As hh itself is very 
small and also because it is very precise, since all of the 
energies involved have tiny statistical errors. We take a 
50% error on b quark mistuning but a 200% error from 
c quark mistuning to allow for the fact that we may be 
underestimating the slope with c quark mass because of 
discretisation errors in the HISQ action for c on the very 
coarse lattices. 



Table VII gives tuned values for A_B c ,hh on each en- 
semble, along with three errors; that from statistics, from 
tuning and from the uncertainty in the lattice spacing. 
These latter two errors dominate. As before, variation in 
the value of the lattice spacing means that masses must 
be retuned. Here this has the effect of producing a net 
change equal to the naive lattice spacing error. Figure [7] 
plots these results against the square of the lattice spac- 
ing, after reconstructing the B c mass by adding back in 
(M b i tPhys + M„ o , phys )/2 = 6.2175 GeV. 

Clear lattice spacing dependence is visible in Figure [7] 
but there is no sign of sea quark mass dependence. As 
for the case of Ab s we fit the results for Ab C) wi as a func- 
tion of lattice spacing and sea quark mass to extract a 
physical result. The fit form is essentially the same as for 
Ab 3 - However, because As c ,hh is such a small quantity 
it cannot set the scale for the discretisation and sea quark 
mass effects. So, instead of allowing a function of a, Sxi 
and Sx s to multiply Ag^j,/, we add such a function with 
a multiplicative factor of 0.4 GeV, representing a typical 
scale for QCD binding energies. We include more terms 
for discretisation errors than in the B s case and set their 
scale by m c , rather than A, because in this case discreti- 
sation errors will come largely from the HISQ action for 
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FIG. 7: Results for the mass of the B c meson tuned to the 
correct c and 6 quark masses on each ensemble obtained from 
the 'hh method' and plotted against the square of the lattice 
spacing. The errors on the points include statistical and tun- 
ing errors and the uncertainty in the lattice spacing. The dark 
shaded band is our physical result, allowing for sea quark mass 
and lattice spacing dependence as described in the text, and 
including electromagnetic effects. The width of the lighter 
shaded band reflects our full error as given in the error bud- 
get, Table |V| The black star is the experimental result [20] . 
offset from a = for clarity. 



the c quark. 



(21) 



&B c Mh(a,Sxi,Sx s ) = A BcM h, P h ys + 

4 

0.4 [ 2J Cj(m c a) 2j ( 1 +c jb 5x m + c jbb (Sx m ) 2 ) 

+2bi6xi( 1 +ci(m c a) 2 + c i; (m c a) 4 ) 
+2b s 5x s ( 1 +c s (m c a) 2 + c ss (m c a) 4 ) 
+4b u (Sxi) 2 + 2b ls 6xi5x s + b ss {6x s ) 2 ]. 

We take the prior on AB c ,hh, P hys to be 0.05(5). The priors 
on the sea quark mass dependence, bi and b s , are taken 
to be 0.00(7), and on the parameters corresponding to 
the quadratic sea mass dependence 0.000(13), as before. 
We take the scale of the physical et-dependence to be the 
scale of m c ~ 1 GeV, as discussed above. The coefficients 
of the quadratic a dependence, c±, c\ and c s , should be 
0(a s ) so we take priors of 0.0(5). For c 2 , and the am b 
dependence of the discretisation errors, Cj b and Cj bb , we 
take a very conservative prior of 0(1). 

The result for A Bci hh,pbys is 0.0616(42) GeV, giving a 
B c mass of 6.279(4) GeV. The fit is robust under changes 
of the prior values. For example, we tried the following 
changes: 

• taking the multiplier of a- and m sea -dependence to 
be 0.8 instead of 0.4; 

• taking the priors for a-dependence to have width 2 
instead of 1; 
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Set 


ami 


om c 


aM Va 


aE Bc 


aAi yp 


arris 


aMo s 


1 


3.4 


0.85 


2.27031(16) 


1.34917(27) 


0.0324(2) 


0.066 


1.5138(7) 






0.85 








0.080 


1.5295(10) 




3.4 


0.88 


2.32148(14) 


1.37456(27) 


0.0325(2) 


0.066 


1.5441(14) 






0.88 








0.080 


1.5587(10) 




3.6 


0.85 




1.35415(29) 


0.0309(3) 








3.6 


0.88 




1.37593(29) 


0.0311(3) 








3 4(c, ^ 1) 


0.85 




1.34987(27) 


0.0323(2) 






2 


3.4 


0.85 


2.26964(17) 


1.34834(34) 


0.0326(3) 


0.066 


1.5140(8) 


3 


2.8 


0.65 


1.84949(11) 


1.11727(13) 


0.0268(2) 


0.0537 


1.2260(5) 


4 


2.8 


0.66 


1.87142(12) 


1.12783(25) 


0.0271(4) 


0.05465 


1.2406(5) 


5 


1.95 


0.43 


1.31691(7) 


0.81861(12) 


0.0210(2) 


0.0366 


0.8709(3) 



TABLE VI: Results for energies and masses needed for the determination of the mass of the B c meson. All the results are in 
units of the lattice spacing. For each valence b quark mass the kinetic mass and energy for the spin-average of T and r\h states is 
given in Table [TTT] as are the r) s meson masses for each s quark mass and the corresponding B s meson energies. In column 4 we 
give the rj c meson mass corresponding to each value of am c , and in the final column we give the corresponding D 3 meson mass. 
Note that the values for M Vl . and Md b are different from those reported in [12] because here we are using a nonperturbatively 
determined Naik coefficien t as discussed in the text. In column 5, we give the fitted energy of the B c meson (i.e. Eq from fits to 
the form given in equation 12 1. In column 6 we give the hyperfine splitting, A^ vp = E(B*) - E(B C ), discussed in section 
This column is largely from 23 but includes some additional values on set 1 that we use for studying systematic errors. 



HID 



Set A Bc , hh (GeV) -A Be , hs (GeV) 
1.034(2) (16) (4) 
1.035(2)(28)(4) 
1.044(3) (9) (4) 
1.046(3)(10)(5) 



1 0.0980(4) (12) (8) 

2 0.0974(4)(35)(9) 

3 0.0782(2) (5) (6) 

4 0.0788(3) (9) (7) 



5 0.0652(3)(14)(5) 1.054(2)(2)(4) 

TABLE VII: Results for Aa c , the mass difference between 
the B c meson and a particular reference mass, on different 
ensembles after tuning to the correct valence b, c and, where 
appropriate, s quark masses. The 3 errors listed are statis- 
tical, tuning and from the uncertainty in the lattice spacing. 
Column 2 gives results from the hh method and column 3 
from the hs method, as described in the text. 



• taking the priors on sea quark mass dependence to 
be 0.5 rather than 0.2. 

None of these changed the result by more than 1 McV. 

Our result is for a world without electromagnetism or 
charm quarks in the sea. The effects of electromagnetism 
on the B c can bc estimated from a potential model in the 
same way as we have done for bottomonium and charmo- 
nium. The quark and antiquark in the B c have the same 
sign of electric charge, however, and so the effect now is 
repulsive rather than attractive. We estimate that the 
effect of switching on electromagnetism is to shift our B c 
mass upwards by 2 MeV. The effects of c quarks in the 
sea can b e estimated following the discussion in subsec- 
tion 



II C as approximately 1 /60 of the hyperfine splitting 



in the B c system, or 1 MeV. This effect is attractive and 
so will counteract the effect of electromagnetism. We 
take the net shift in the B c mass as 1(1) MeV, moving 
our result to 6.280(4) GeV. This is the value given by the 
dark shaded band in Figure [7| 

The complete error budget is given in Table [V] Here 



we have split up the 4 MeV error from the fit into its 
component parts as discussed in the B s case and added 
sources of systematic error. Errors from missing relativis- 
tic corrections to the NRQCD action are similar to those 
for the B s case. The leading missing spin-independent 
terms in the NRQCD action are 0(a s v^) and 0(v%). We 
must estimate the effect of these terms on both M fe r and 
M b c . vl is about half the size in the B c compared to bot- 
tomonium. This means that there is some cancellation of 
the otsVfr errors in AB a ,hh, since we estimated a s v^ errors 
previously as a s v^(500 MeV). Independent confirmation 
that these terms have a small net effect comes from the 
calculations that we have done here for the case where C\ , 
C5 and cq take the values that include the 0(a s ) radiative 
corrections, see Table VI and Figures [5] and [6] There will 



be little cancellation of the v% errors, however, since v% is 
much smaller in the B c than in bottomonium. We take 
the systematic error from spin-independent terms then to 
be 4 MeV for a s vf (i.e. half that for B s ) and 5 MeV for 
v%, added in quadrature to give 6 MeV. The leading spin- 
dependent error is from missing radiative corrections to 
the a.B term. This affects the B c only because of the 
spin-averaging of M b ^. As for the B s we take 3/4 of a 
possible 10% error in the hyperfine splitting, estimated 
in [52] at 53 MeV, i.e. 4 MeV. 

Errors from the uncertainty in M bb and M rjc do not 
affect As c ,h/i significantly but Mg a inherits an error of 
half their uncertainty when it is reconstructed from Ab c 
and those masses. These two uncertainties are partly cor- 
related, because they both contain estimates of electro- 
magnetic and annihilation effects in the two very similar 
charmonium and bottomonium systems. The systematic 
errors from electromagnetism and c quarks in the sea 
for the B c are also correlated with the errors for these 
effects from charmonium and bottomonium. These er- 
rors are marked with a star in Table [V] and we separate 



13 




0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 
M lls 2 (GeV 2 ) 




(GeV) 



FIG. 8: Results from set 1 for the mass of the B c meson 
(specifically the difference between that mass and the masses 
of the B„ and D s mesons) as a function of the square of the 
r/ s meson mass, acting as a proxy for the strange quark mass. 
The errors are statistical only, since lattice spacing errors af- 
fect all the points together. The lines are fits to the results 
allowing linear terms in M% s , M hh - and M Va . Here the lines 
join points for fixed b and c quark masses. See Figure [9] for 
the equivalent as a function of M bb and Figure [10] for the 
equivalent as a function of M Vc . 



FIG. 9: Results from set 1 for the mass of the B c meson 
(specifically the difference between that mass and the masses 
of the B s and D s mesons) as a function of the spin-average 
mass of the T and r/b, M bb , acting as a proxy for the b quark 
mass. The errors are statistical only, since lattice spacing 
errors affect all the points together. The lines are fits to the 
results allowing linear terms in M^ s , M bb and M Vc . Here the 



lines join points for fixed s and c quark masses. See Figure [8] 



for the equivalent as a function of Mi and Figure 
equivalent as a function of M Vc . 



10 



for the 



out the different components and add them linearly with 
appropriate signs before squaring and accumulating into 
the total error. The error from these three components 
is then increased by their correlation from 3 MeV to 4.6 
MeV. We estimate finite volume errors to be negligible for 
the B c . Our total error, from adding statistical and sys- 
tematic errors in quadrature, is 9.5 MeV, giving a mass 
for the B c meson from the hh method of 6.280(10) GeV. 
The total error is plotted as the more lightly shaded band 
in Figure [7} 

The hs method has different systematic errors from the 
hh method and so provides a good cross-check. The raw 



results needed for this method are given in Table VI and 
in Figures [HJ [9] and [lO] we show results for the quantity 
(because Ag c j is is negative) defined by: 



A Bo , hs = {E Bc - (E Bs +M Ds )) 



latt ' 



(22) 



In the figures —Ab c ms results from set 1 are plotted 
against the different quark masses involved in the cal- 
culation, with M^, M b ^ and M, lc acting as proxies for 
the s, b and c quark masses respectively (we have two 
different values for the masses of each quark). We see 
that there is fairly strong dependence on the s quark 
mass but very little on the b quark mass or the c quark 
mass. The slope against is 0.41, which agrees well 
with that expected from experiment if we substitute light 
quarks for s quarks in the formula for Ab c: /, 8 above. The 
slope against is 0.005 and against M Vc is 0.07. These 
are only in very rough agreement with the expectations 




FIG. 10: Results from set 1 for the mass of the B c meson 
(specifically the difference between that mass and the masses 
of the B 3 and D s mesons) as a function of the mass of the 
rj c acting as a proxy for the c quark mass. The errors are 
statistical only, since lattice spacing errors affect all the points 
together. The lines are fits to the results allowing linear terms 
in M^ 3 , M bb and M Vc . Here the lines join points for fixed s 
and b quark masses. See Figure [8] for the equivalent as a 
function of 
of M hb . 



and Figure [9] for the equivalent as a function 
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FIG. 11: Results for the mass of the B c meson tuned to the 
correct c and b quark masses on each ensemble obtained from 
the 'hs method' and plotted against the square of the lat- 
tice spacing. The errors on the points include statistical and 
tuning errors and the uncertainty in the lattice spacing. The 
darker shaded band is our physical result, allowing for sea 
quark mass and lattice spacing dependence as described in 
the text, and including electromagnetic effects. The width 
of the lighter shaded band reflects our full error as given in 
the error budget, Table [V] The black star is the experimental 
result [20], offset from a = for clarity. 



of 0.02 and 0.04 respectively from comparing results over 
the much larger experimental range from b to c |20j . 



As before we can use these results to compensate for 
mistuning of the quark masses. Again we take a 50% er- 
ror on tuning shifts for b and s but a 200% errors on those 
for c to allow for discretisation errors in our estimates of 
those effects. The tuned results are given in Table |VII| 
All of the statistical, tuning and lattice spacing errors 
are larger than those of the hh method. The tuning er- 
ror dominates on the very coarse and coarse lattices, but 
on the fine ensemble it is comparable to the other er- 
rors. The lattice spacing error is reduced by a factor of 
two over the naive error by the retuning required when 
the lattice spacing changes. Lattice spacing dependence 
is small but visible in these results; no sea quark mass 
dependence is evident. 



As for the B s and the B c hh method we fit the re- 
sults for As ci / ls to a functional form allowing for lattice 
spacing dependence, including NRQCD effects, and sea 
quark mass dependence. The functional form is the same 
as that used for the hh method, except that now we can 
include these dependences as a multiplicative factor since 



Ab c ,/is is not unusually small. We use: 

^■B c ,hs(a,Sxi,Sx s ) = A B<; ,h S) phys [l + (23) 

4 

^c.,(TO c a) 2j ( 1 +c. 

3=1 

+2biSxi( 1 +ci(m c a) 2 + c u (m c a) 4 ) 

+2b s 5x s ( 1 +c s (m c a) 2 + c ss (m c a) 4 ) 

+Ab u (8xi) 2 + 2bi s 5x l Sx s + b ss (5x s ) 2 ]. 

We take the same prior values as for the hh method ex- 
cept for the prior for Ab hs phys which we take to be 
-1.0(2). 

We obtain the result A BeihSiPhys = -1.054(17) GeV. 
This changes by less than 1 MeV on doubling the prior 
width for the lattice spacing dependence or the sea quark 
mass dependence. To reconstruct the B c mass from this 
we need to add the appropriate values for the B s and D s 
masses in a world without electromagnetism or c quarks 
in the sea. As discussed earlier, electromagnetism has 
negligible effect on the B s mass. The D 8 mass is lower 
by 1.3(7) MeV, however, in a world without electro- 
magnetism, from the phenomenological formula in equa- 
This gives a total for the appropriate sum of 



tion 17 

Md 3 



of 7.334 GeV. Figure 11 shows our tuned 
results for ms a as a function of lattice spacing. The dark 
shaded band is the result from the fit just described in- 
cluding the error obtained from it. We have shifted the 
B c mass obtained upwards by 1(1) MeV, to a central 
value of 6.281 GeV, as described earlier to allow for elec- 
tromagnetic and charm-in-the sea effects that are not in- 
cluded in our calculation. The lighter shaded band gives 
the total error, from the error budget of Table [V| the sys- 
tematic error components of which we will now discuss. 
The first three entries in the Table are the split of the 
17 MeV error obtained from the fit among its different 
components. 

Errors from missing relativistic corrections to the 
NRQCD action affect the B c energy and the reference 
B s energy. The leading missing spin-independent terms 
in the NRQCD action are 0(a s v%) and 0{vf). We ear- 
lier estimated the shift from O{ot s vf) terms on M.g at 
15 MeV. For the B c we expect a systematic error of 
about half this value, so we take 7.5 MeV, since v\ is 
roughly half as big. For B s a s vf terms have very little 
effect and neither meson will be sensitive to terms. 
Spin-dependent NRQCD errors come chiefly from miss- 
ing radiative corrections to the a ■ B term, but there will 
be cancellation here between the B c and the B s since 
both mesons will respond in the same way to a change 
in C4. We therefore take a systematic error of 1 MeV 
rounding up the difference between the 4 MeV system- 
atic previously allowed for this for the B c and the 3.5 
MeV systematic for the B s . 

Systematic errors from uncertainties in the physical 
values of M Vb , M Vc and M b i which affect the quark mass 
tuning can be estimated from the dependence of Ag^^s 
on these quantities discussed earlier. The error from the 
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uncertainty in M rja is sizeable at 2.3 MeV; the others 
are very small. We must also allow for systematic errors 
from uncertainties from electromagnetism and charm-in- 
the-sea for the reference masses of the D s (0.7 MeV is half 
the shift applied in that case) and the B s (negligible) as 
well as for the B c itself (1 MeV as discussed above). 

This gives a total error of 19 MeV, dominated by the 
statistical and tuning errors of the raw data. Our fi- 
nal result for the B c mass from the hs method is then 
6.281(19) GeV. This is plotted as the more lightly shaded 
band in Figure [TTJ The agreement between the hs and 
hh methods is very good, although their systematic and 
statistical errors are very different, with the hh method 
being significantly more accurate. The agreement is in 
fact not surprising when we consider that the B s mass 
determined in section [ill A| agreed well with experiment. 
The B c hs method replaces M b ^ with Mb 3 and M Vc with 
Md 3 so if the B s and D s masses are known to agree 
with experiment given masses tuned from M.t and M„ 
then the B c from hs and hh will agree. However, the 
fact that they were derived completely independently is 
a good consistency check of the method and of our error 
estimates. 



B mass 



0.12 
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FIG. 12: Results for the mass difference between the B a and 
B mesons as a function of Ml — M% acting as a proxy for 
m s — mi. Open squares are results on very coarse ensembles, 
triangles on coarse ensembles and the cross is from the fine 
ensemble. The green band represents the result of our fit, 
adjusted for electromagnetic effects as described in the text. 
The black star is the experimental result for Mb b — (M B o + 
M B ±)/2. 



The correlators for the B meson are noisier than those 
for the B s , as will be clear from the discussion of the sig- 
nal/noise earlier. This means that the B meson mass is 
the least well determined of all our masses. The best 
way then to pin down the B mass is to consider the 
mass difference between the B s and the B. NRQCD 
systematic errors should entirely cancel in such a differ- 
ence. However, because the difference is a small number 
we have a fairly sizeable statistical error even when we 
fit both correlators together, as described earlier, and ex- 
tract Eb, — Eb directly from the fit. 



Tabic VIII gives values for the energy difference be- 
tween B s and B extracted from our fits on each ensemble. 
Statistical errors are 10-15% of the splitting. However, 
this amounts to less than 10 MeV in terms of the absolute 
mass, so still provides a good test against experiment for 
M B . 

We plot the results for M Bs —M B (= a~ 1 (aE Bs —aE B )) 
against M? — M 2 which is a useful physical proxy for 
m s — mi in Figure 12 We expect this mass difference to 
be largely linear in m s — mi and our results are consistent 
with this. Given the statistical errors, we fit a relatively 
simple form to this difference: 

(M Ba - M B )(a,5xi,5x s ) = 

3 

5^ ajfj{M n „ M„)(l + C]1 (Aa) 2 + c j2 {Aa) 4 ) 
i=i 

+f l (M 11sl M J! ){2b l 8xi{l + ci{ka) 2 )) 

+2b s Sx s {l + c s {Aa) 2 ). (24) 



Here the functions fj are simple ones that respect the 
fact that Mb, ~ Mb vanishes by definition when the 
light quark mass is equal to the strange quark mass. So 
fx = M 2 s - Ml f 2 = Ml - M* and / 3 = M* log M„; - 
M 2 log M 2 . We allow these terms to have lattice spacing 
dependence with a scale set by A = 400 MeV. We also 
allow sea quark mass dependence in the terms multiply- 
ing j\. The coefficients aj are given priors of 0.0(5) (we 
expect a slope of 0.2 if the dependence on M 2 s —M 2 were 
purely linear). For the Cj% a 4 -dependence coefficients we 
take 0.0(1.0) and for cji we take 0.0(5) since a 2 terms 
should be suppressed by an additional power of a s . For 
the sea quark mass dependent coefficients bi and b s we 
take priors of 0.00(7) as discussed earlier. 

The physical value for Mb s — Mb is then the value 
at M Vs = 0.6858 GeV and M w o = 0.135 GeV, for a = 
Sxi = 8x s = 0. We obtain 0.073(14) GeV. This value 
is to be compared to the experimental mass difference 
between the B s meson and the average of the charged 
and neutral B mesons which we denote Bi (thus aver- 
aging the u and d quark masses). However it has to be 
adjusted for electromagnetic effects not included in our 
lattice QCD calculation. Following the discussion in sec- 
tion III A we see that electromagnetic effects in the B s 
and Bd mesons are very small. For the B u however, be- 
cause it is charged, the shift is more substantial at 2.2 
MeV Adding in electromagnetism then shifts our 

M Bs - M Bl splitting down by 1 MeV. The result 72(14) 
MeV is shown as the shaded green band in Figure [l2j It 
is in reasonable agreement with the experimental result 
of 87 MeV [D]. 
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Set 




am 3 


ami aMir 


aE{B 3 ) - aE(Bi) 




1 

2 


3.4 

Q 1 

3.4 


0.066 
0.080 
0.066 


0.0132 0.2408(6) 
0.0132 

0.0264 0.3348(6) 


0.0553(62) 
U.Uooo^oz ) 
0.0369(29) 


0.0318(78) 
0.0374(55) 


3 
4 


2.8 
2.8 


0.0537 0.0067 0.1567(4) 
0.05465 0.01365 0.2222(5) 


0.0446(46) 
0.0336(41) 


0.0306(54) 
0.0245(68) 


5 


1.95 


0.0366 


0.00705 0.1377(4) 


0.0245(30) 


0.0177(35) 



TABLE VIII: Results for energies and masses needed for the determination of the mass of the B meson. Column 6 gives the 
energy difference between B s and Bi mesons, for different valence b, s and I quark masses given in columns 2, 3 and 4. Column 
5 gives the corresponding n meson mass, taken from Column 7 gives the hyperfine splitting for the Bi meson, from [22], 
A hy P = E{B!) _ E{m y 



Our final result for M Bl = (M B ± + M B o)/2 is 5.363 - 
0.072 = 5.291(11)(14) GeV. The first error comes from 
the mass of the B s and is discussed in subsection |III A[ 
the second comes from the mass difference between the 
B s and Bi. We do not expect any significant additional 
systematic errors from NRQCD, beyond those that the 
Bi inherits from the B s in this method, because those 
errors should cancel in M Bg — M Bl . The error budget 
for M B[ is then as given for M Bs in Table [V] with the 
additional 14 MeV given above. Our result for M Bl of 
5.291(18) GeV can be compared to the experimental re- 
sult of 5.2794(3) GeV [2UJ. 



D. B* - B splittings 

By projecting out the vector states at the source and 
sink we can measure the correlator for the B* , B*, and 
B*. As they come from exactly the same configurations 
and valence HISQ propagators as the corresponding pseu- 
doscalar states, they are highly correlated. In this case 
we can do simultaneous fits to both the pseudoscalar and 
vector meson propagators, and extract a value for the 
M B * — M B splittings. 

This hyperfine splitting is generated by the a ■ B term 
in the NRQCD action, equation [6] This is 0(v^) in the 
relativistic power counting for heavyonium and C(A/Mf,) 
in heavy-light power counting. In our action we take the 
coefficient of this term, C4 to be 1, but it will have radia- 
tive corrections when matched through 0(a s ) with full 
QCD. We are also missing higher dimension operators 
that correct for discretisation errors and add relativistic 
corrections. For B systems, which are relatively large 
with very slow-moving b quarks, we do not expect these 
latter effects to be as important as the issue of the deter- 
mination of C4 beyond tree level. The heavy-light hyper- 
fine splitting generated by the a ■ B term is proportional 
to C4 and so uncertainty in C4 leads directly to an 0(a s ) 
i.e 25%, uncertainty in the splitting which decreases only 
slowly on finer lattices. Thus to determine this splitting 
accurately we need a determination of C4. 

Since we use exactly the same NRQCD action for all 
our calculations, however, we can effectively determine C4 
by comparing one set of heavy-light hyperfine splittings 



to experiment and then predicting the others. Equiva- 
lently we can take ratios of hyperfine splittings in which 
the normalisation factor, C4, cancels. This is what we 
did in [22]. By using the B* — B s mass difference, which 
is 46.1(1.5) MeV from experiment [20], we showed that 
this splitting does not depend on the mass of the lighter 
quark even for as heavy a quark as the charm quark, and 
we were able to predict a B*—B c splitting of 53(7) MeV. 

We will not discuss that analysis further here, but we 
give the table of results of the hyperfine splittings for 
completeness in Tables inj VlJ VIII| They include some 
additional values over those in [22 for the purposes of 
further testing systematic errors. Figure [13] shows such a 
test in a plot of the mass difference between B* and B s 
as a function of M,-r for two different s quark masses (the 
first 4 rows of entries in Table III ) . Dependence on the 
b quark mass is visible, but dependence on the s quark 
mass is very small. Results are also shown for B* — B c 
for the same b quark masses and they show a parallel 
slope. In addition we show a result for the the NRQCD 
action with a coefficients different from 1 and no change 
is seen. 

Here we are interested in analysing the systematic error 
in the B, B s and B c meson masses from any uncertainty 
in C4. We do this by comparing our B* — B s splitting 
to experiment and interpreting any mismatch as a signal 
for C4 ^ 1. We use the B* — B s because this is the most 
accurately determined splitting from our analysis that is 
also known experimentally. This method can be used as 
a nonperturbative determination of C4, and we used this 
previously to bound the errors on our prediction of the 
hyperfine splitting in bottomonium based on earlier B s 
and B hyperfine splitting results [IT]. Figure 14 shows 
our results as a function of lattice spacing. We have ad- 
justed them for mistunings of the b quark mass according 
to the results in Figure [13] but the corresponding shifts 
are small, and less than the statistical errors in all cases. 
The results show little sign of any lattice spacing depen- 
dence or sea quark mass dependence and we see that a 
value of C4 of w 1.1 would give agreement with experi- 
ment for all the values. We therefore estimate that the 
correct value of C4 for this NRQCD is 1.1(1) and that we 
make a systematic error of about 10% in the heavy-light 
hyperfine splitting by using C4 = 1. This produces a sys- 
tematic error on the B c and B s meson masses discussed 
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FIG. 13: Results for the mass difference between the £?J and 
B s mesons (open squares) and between B* and B c mesons 
(crosses) as a function of the spin-average of the T and r\h 
meson masses used to tune the b quark mass. BI" — B s results 
include two values of the s quark mass, and B* — B c results 
include two different values of the c quark mass as well as 
results for ci, c 5 and eg coefficients differing from 1. 



FIG. 15: Results for the mass difference between the scalar 
B* and B s mesons as a function of lattice spacing. The 
squares show our results (with multiple b and s quark masses 
on the very coarse ensemble, set 1) and the shaded band the 
physical result from the fit described in the text. This result 
does not include any adjustment or error for the fact that the 
scalar is not a gold-plated meson. 
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FIG. 14: Results for the mass difference between the BI" and 
B B mesons as a function of lattice spacing. Results have been 
corrected for mistuning of the b quark mass and the errors 
include statistics, mistuning and lattice spacing uncertainties. 
The black star and block dashed lines give the experimental 
result [20]. 



in earlier subsections and included in Table |V] 

Note that the behaviour of the hyperfine splitting in 
bottomonium is quite different from that of the B s , being 
strongly dependent on the lattice spacing [TT]. However, 
it is the same operators in the NRQCD action, with the 
same coefficients, that control the fine structure in both 
systems. The matrix elements of the operators can be- 
have quite differently, and bottomonium is expected to 
be a lot more sensitive to the lattice spacing than the 



B s . This means that the B s is a good system from which 
to determine C4 because it is really only sensitive to that 
coefficient. 



E. Scalar and Axial vector meson masses 

When generating the NRQCD propagators we choose 
Dirac structures T to explicitly project out pseudoscalar 
and vector mesons. Parity partners of both of these con- 
tribute to their correlators, as shown in equation [l2j and 
must be included in the fit. The parity partner state of 
the pseudoscalar is a scalar meson and the vector meson 
has as its parity partner a axial-vector state. So, by care- 
fully fitting the correlators of the 0~ and 1~ states we 
get also the spectra of the + and 1 + states for free. 

In fact our fit results return directly the mass difference 
between the ground state in the oscillating channel and 
the ground state in the non-oscillating channel i.e the 
+ — and 1 + — 1~ mass differences. We report these 

(i.e. 



for the B* s0 and B sl 



results in Table III for the B, 
mesons) . For B\ and B c our errors are too large on some 
fits to give a full picture across all ensembles. 

The results for the scalar-pseudoscalar mass difference 
are shown in Figure 15 There is no signal for any sys- 



tematic dependence on the b or s quark mass, or on the 
lattice spacing. In deriving a physical result we allow for 
both physical and unphysical lattice spacing dependence 
as described for the B s mass in subsection III A[ as well 
as sea quark mass dependence. We use the same fit form 
as for the B s mass, given in equation |12| The priors are 
taken to be the same except that we take the prior on 
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FIG. 16: Results for the mass difference between the axial 
vector B s i and B* mesons as a function of lattice spacing. 
The squares show our results (with multiple b and s quark 
masses on the very coarse ensemble, set 1) and the shaded 
band the physical result from the fit described in the text. 
This result does not include any adjustment or error for the 
fact that the axial vector is not a gold-plated meson. 



the physical value of the scalar-pseudoscalar mass split- 
ting to be 0.4(2). We also allow for more sea quark mass 
dependence than in that case, because the scalar meson 
is not gold-plated (this will be discussed further below). 
We therefore do not take the factor three suppression of 
sea quark mass effects in this case, so the prior on the 
sea quark mass dependent terms is simply 0.0(2) for the 
linear terms and 0.00(4) on the quadratic terms. The 
physical result we obtain is 0.385(16) GeV and this is 
plotted as the shaded band on Figure \TE\ 

Exactly the same procedure is followed for the axial 
vector - vector splitting. The results are plotted in Fig- 
ure [16] From the same fit as that described above we 
obtain the physical result for the mass difference between 
the axial vector and vector of 0.391(15) GeV, plotted on 
the Figure as a shaded band. 

Since we have calculated the mass differences between 
the scalar and axial vector B s mesons and the corre- 
sponding pseudoscalar and vector B s mesons we expect 
only very small systematic errors coming from NRQCD. 
Because the b quark is very nonrelativistic in these sys- 
tems, as discussed in subsection |III A| the errors from 
missing higher order relativistic corrections in NRQCD 
are very small. They will be reduced further here by 
cancellation in the mass difference. The main source of 
systematic error from NRQCD will come from radiative 
corrections to spin-dependent terms in the NRQCD ac- 
tion. In subsection |III D| we showed that these systematic 
errors are not large, at least for the a ■ B term. There the 
errors amounted to 10% of the hyperfmc splitting, around 
5 MeV, split between the vector and pseudoscalar states. 
Assuming a similar error for other spin-dependent terms 



which would affect p-wave states, we take a systematic 
error of 5 MeV from NRQCD in the mass differences. 

A potentially much larger source of systematic er- 
ror is the fact that the scalar and axial vector mesons 
have strong decay modes, i.e. they are not 'gold-plated'. 
This will be discussed further in the next section. The 
strongest decay mode, if kinematically allowed, will be to 
BK (for the B* ) or B*K (for the B sl ). If the masses are 
such that the mesons are below threshold for this decay 
mode, there will still in principle be coupling between 
the meson and this virtual decay channel which can shift 
the meson mass. There is in addition a Zweig-suppressed 
decay mode to B s tt/B*t: which will be kinematically pos- 
sible. 

On the lattice the coupling between single and mul- 
tiparticle states is distorted by the fact that u/d quark 
masses are heavier than their physical values and the vol- 
ume of the lattice is relatively small. The fact that u/d 
masses are unrealistic means that decay thresholds are 
higher than in the real world. In principle sensitivity to 
decay thresholds would be seen in the results as sea quark 
mass dependence, but that may not become visible until 
much closer to real world u/d mass values. The finite vol- 
ume of the lattice restricts the decay momenta that real 
or virtual multiparticle states can have. A lattice anal- 
ysis on multiple volumes allows single and multiparticle 
states to be separated. In practice [3T] it seems that bi- 
linear operators of the kind that we have used here have 
very small overlap with multiparticle states. So, although 
in principle there may be a multiparticle state (such as 
B s n) at a lower mass value than the B s q it is very hard 
to pick it out of a lattice QCD calculation without ex- 
plicitly using multiparticle operators, which we have not 
done. 

A simple model to analyse the effect of multiparticle 
states is one in which point-like meson states are coupled 
together via a perturbation which is a simple point-like 
vertex. We can then calculate the shift on the single 
particle energy from this coupling by integrating over the 
momenta of the decay products in the initial particle rest 
frame. For example, for B* coupling to BK: 



A, 



d 3 p 



1 



(2n)3M B , Q -(E B +E K ) 



(25) 



A represents an ultraviolet cut-off required for this model 
to make sense. A w 500 MeV and g 2 w (0.5/A). If our 
calculation is correct that the B* is close to, but below, 
threshold then we can treat the B and K as nonrelativis- 
tic and, dropping the B kinetic term, 



~2tt 



p 2 dp 



p 2 /(2M_ 



K J 



AM ' 



(26) 



where AM = Mb + Mk — Mb; (unperturbed values). 
The AM-dependent piece of the mass shift is then given 
by: 



A, 



g 2 y/AM2M K 2M K 



(27) 
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Numerically this gives a shift downwards of a few tens 
of MeV for AM values of a few tens of MeV. From this 
we conclude that a reasonable systematic error for the 
absence of coupling to strong decay channels is 25 MeV 
(which we take to be a symmetric error). 

This then gives the following mass differences: 



M B * o -M Bs = 0.385(16) (5) (25) 
M Bsl -M B , = 0.391 (15) (5) (25) 



(28) 



where the first error is statistics/fitting, the second is the 
NRQCD systematic error and the third is the error from 
not including coupling to strong decay channels. 



IV. DISCUSSION 

Figure [4] shows that our result for the mass of the B s 
meson agrees well with experiment with total errors of 
11 MeV (0.2%). The errors are dominated by statistical 
errors and systematic errors from NRQCD, both of which 
are being improved in work underway. 

As discussed earlier, because we fix the b quark and 
s quark masses from other mesons, the B s mass deter- 
mination is completely free from any parameter tuning. 
An alternative for the b quark mass, adopted by some 
other lattice QCD calculations is to fix the b quark mass 
from the B s meson mass itself. However, it is still possi- 
ble then to determine = M Bs — M b r/2, as a test of 
the b quark systematic errors. The only other full lattice 
QCD calculation of this quantity is from the Fermilab 
Lattice/MILC collaboration using the Fermilab formal- 
ism for the b quark [25) . They determine in fact the 
quantity A^- = Mg- — M bb /2 where B s is the spin av- 
erage of the B s and the B* masses. This quantity has 
reduced systematic errors from the spin-dependent terms 
in the action, in the same way that the use of M bb reduces 
these systematic errors for the bottomonium system. The 
Fermilab Lattice/MILC collaboration obtain the value 
1359 ± 304+^ MeV for 2A^ with a partial error bud- 
get [35] . We can also determine A^- in exactly the same 
way as we determined A^. We obtain 0.671(7) GeV for 
the physical result from our calculation. This becomes 
0.675(11) GeV when corrected for electromagnetic, an- 
nihilation and charm-in-the-sca effects in bb and with a 
full error budget (essentially the same as in Table [v] but 
with a reduced error for NRQCD systematics in the B s ). 
The experimental result is 0.6817(11) GeV [2D] . 
Figure 
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shows the results for A Bs and A-g- from this 
paper andfrom the Fermilab Lattice/MILC collabora- 
tion compared to experiment. Both results agree with 
experiment but we are able to provide a 2% test of these 
mass differences, which is a nontrivial test of QCD. 

Back in 2004 we predicted the mass of the B c ahead 
of the CDF experimental discovery in a collaboration 
with the Fermilab Lattice collaboration 26J. We used 
NRQCD for the b quarks, as here, but the Fermilab for- 
malism for the c quarks, and the asqtad formalism for the 
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FIG. 17: Results for the mass differences Ab s and Ag- be- 
tween the B s and the spin average of B 3 and B* respectively 
and the spin average of the T and r/b (see text). The top 
result is from the Fermilab Lattice/MILC collaboration for 
Ag- [25], the middle two results are from this paper, and the 
lower two points, and shaded vertical lines, are from experi- 
ment. 
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FIG. 18: Results for the mass of the B c obtained from the 
hh and hs methods. The top two results (hh above hs) are 
the 2004 HPQCD/Fermilab Lattice calculation [2H] and the 
bottom two results (again hh above hs) are from this paper. 
The middle point, and the shaded vertical line, is the current 
experimental value. 



20 



700 



(,00 



300 



200 



■ 






- 


• 
• 

• 


• 

• • 




. : . : 

: 




• 


$ 


• 


• 










D 




b xA x b - 



FIG. 19: Results for the mass difference between + , 1 + and 
2 + D, D s , B and B s states and the spin-average of the cor- 
responding 0~ and 1~ states. The solid circles are from ex- 
periment [20] and are divided, where results exist, into the 
expected + /l + and l + /2 + doublets. For the B and B s only 
the upper l + /2 + doublet has been seen. The open circles 
for the B s are the lattice QCD results given here. We have 
not included an error for the coupling of these states to de- 
cay channels. The rightmost two sets of points are the P- 
wave charmonium and bottomonium states for comparison. 
The black lines show the D M K and B M K thresholds for the 
Zweig-allowed decay of the + and 1 + D s and B s states. 



s quarks in the hs method. As a result, we had larger 
statistical and systematic errors than we have here, par- 
ticularly for the hs method. 



Figure 18 shows the comparison between our old re- 
sults and the new ones given here, as well as the current 
experimental value. The improvements in lattice QCD 
calculations since 2004, including the development of the 
HISQ action for c and s, give us a substantial improve- 
ment in errors and consistency between the hh and hs 
methods today. 

Our result for the mass difference between the B s and 
the B[ meson is the first full lattice QCD calculation of 
this quantity. As discussed in subsection III C our result 
is in agreement with experiment, but with substantial 
statistical errors. These will be improved in further work 
which is under way. 

In subsection |III E| we gave results for the masses of 
the + and 1 + B s mesons. The + has not been seen 
experimentally. A 1 + state has been seen but may not 
be the one whose mass we have calculated. 

Figure [19] shows how our results fit into the current ex- 
perimental picture of 'p-wave' charm-light and bottom- 
light mesons. Charmonium and bottomonium p-wave 
mesons are also shown for comparison |20j . For heavy- 
light mesons the p-wave states are expected from heavy 
quark symmetry |27j to appear in two doublets, classified 
according to the J of the light quark which can be either 
1/2 or 3/2, when L—l is combined to s/ = 1/2. The 
jl = 1/2 doublet then separates into a + and 1 + meson 



at non-infinite heavy quark mass, when the heavy-quark 
spin is coupled in. The ji = 3/2 doublet is likewise made 
up of 1 + and 2 + states. This is in contrast to the heavyo- 
nium case where there is a triplet of + , 1 + and 2 + states 
with total quark-antiquark spin of 1, and a single 1 + state 
with total spin 1. The existing experimental results are 
shown as solid points in Figure [19] divided appropriately 
according to the picture above. For D and D s mesons 
both doublets have been seen; for B and B s mesons only 
the ji = 3/2 doublet has been seen (assuming that the 1 + 
state seen is associated with that doublet). For charmo- 
nium all 4 states of the Xc triplet and the h c are known; 
for bottomonium the hy, has not been seen. The exper- 
imental masses are given relative to the spin-average of 
the s-wave states, a pseudoscalar and a vector in every 
case. That removes the overall mass scale of each sys- 
tem from the plot and shows, as is well-known but still 
somewhat surprising, that the orbital excitation energies 
of heavy degrees of freedom in heavyonium are very sim- 
ilar to those of light degrees of freedom in a heavy-light 
system. 

Since mass splittings between the 5=1 states in 
heavyonium and between the members of the ji = 1/2 
or ji = 3/2 doublets in the heavy- light case are caused 
by heavy quark spin effects proportional to the inverse 
of the heavy quark mass we expect to see larger split- 
tings in the c case than in the b case. This is borne out 
in the experimental data for charmonium and bottomo- 
nium and in the comparison of D and B results for the 
ji = 3/2 doublet (although the disagreement between B 
and B s might indicate that the doublet assignment for 
the B 1 + in the Figure is wrong). The splitting between 
ji = 1/2 and 3/2 doublets is a light quark effect that does 
not vanish as rriQ —¥ oo. However the splitting will vary 
with rriQ slightly because of A/toq terms in the effective 
heavy quark action (NRQCD) away from that limit. A 
variation in the splitting of order 100 MeV out of 500 
MeV is then reasonable between D and B. 

Our results are entered on Figure [19] as open circles in 
the B s column. Since we have + and 1 + states we have 
placed them as the j\ = 1/2 doublet. However, it should 
be stressed that we do not know that that assignment for 
the 1 + is correct. In any case the 1 + states from the two 
doublets can mix and we have not allowed for that. 

Given the discussion above, our results fit fairly natu- 
rally into the picture described. As a ji = 1/2 doublet 
they sit below the known ji = 3/2 doublet. They sit 
closer to the ji = 3/2 doublet than for the D s case, but 
this can be a A/rriQ effect as discussed above. The split- 
ting between the two states in the doublet is about one 
third that in the D s case, consistent with this splitting 
being a I/thq effect. 

In the D s case the discovery of the lower + /l + dou- 
blet [2"BH3"U] caused much surprise because the states were 
low compared to model calculations. The states had been 
expected to be above threshold for strong decay to DK 
and D* K respectively and therefore broad (unlike the up- 
per doublet which has to decay in a d-wave) . Instead they 
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are below threshold and so decay to the Zweig-suppressed 
D s tt and Z?*7r channels and are narrow. We have marked 
the DK, D*K, BK and B* K thresholds on Figure [19] 
Like the D* our + B 8 state is also below, but very 
close to, its Zweig-allowed decay threshold. A similar 
situation holds for the 1 + state. This might indicate 
that these states would be narrow. However, there will 
also be effects from coupling to the decay channel that 
are not included in our calculation. On our lattices the 
light sea quark masses are heavier than in the real world 
and hence the K mesons containing a valence s quark 
and a sea light quark would be too heavy to allow the 
B* to decay to BK. B s n decay is allowed but with a 
very restricted phase space compared to the real world. 
Coupling to these channels would in principle show up 
as sea quark mass dependence but would need a bigger 
range of sea quark masses than we have used. In subsec- 
tion [IIIE] we allowed a 25 MeV systematic error for these 
coupled-channel effects, noting that the coupling to BK 
decay will tend to push the mass down. 

Our results are the only ones in full lattice QCD to 
date and with realistic b quarks. There have, however, 
been several recent lattice QCD results for the case of u 
and d sea quarks only and taking b quarks in the static 
limit [3TH34] . The most complete is that of the ETMC 
collaboration [M]. They give a mass difference between 
the B s ji = 1/2 doublet and the spin-average of B s and 
B* of 413(12) MeV, with an estimated additional possible 
systematic error of 20 MeV, including coupling to multi- 
particle states. Using experimental results from charmed 
mesons to estimate l/rrib corrections to the static limit, 
they conclude, as we have done, that the scalar B s state 
is close to the BK threshold. 



V. CONCLUSIONS 

We have given the first accurate result for the B s me- 
son mass from lattice QCD including the effect of u, d 
and s sea quarks, and with a full error budget. We have 
improved significantly on an earlier value for the B c me- 
son mass, achieving smaller errors and better consistency 



between two different methods. The determination of 
both of these masses provides a strong test of our lat- 
tice QCD approach to b physics, because they test the 
consistency of heavyonium and heavy-strange or heavy- 
charm physics from the same heavy quark action. All of 
the QCD parameters used here are tuned from other cal- 
culations so our results are parameter free tests of QCD 
against experiment. 

The mass of the B meson, specifically the difference be- 
tween B s and B meson masses, depends on light quark 
physics since heavy quark effects cancel. Our result 
agrees with experiment but needs higher statistical pre- 
cision for a good test. 

We also discuss scalar and axial vector meson masses 
for the B s . Our results indicate masses below, but close 
to, threshold for Zweig-allowed decay modes. From our 
current calculation, however, it is not possible to include 
effects of coupling to either allowed or suppressed decay 
channels, so significant shifts to our results from these 
effects are possible. 

Further improvement to these results will come with 
improved statistical accuracy in calculations now under- 
way. This will lead also to improved determination of 
decay constants and other B meson matrix elements. 
Confidence in those calculations and the error analysis 
associated with them is strongly bolstered by this analy- 
sis of the associated meson masses. 
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